First principles simulations of the magnetic and structural properties of Iron 
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We have implemented non-coUinear GGA and a generalized Bloch's theorem to simulate uncon- 
mensurate spiral arrangements of spins in a Density Functional Theory code based on localized wave 
functions. We have subsequently performed a thorough study of the different states of bulk Iron. We 
determine the minimal basis set required to obtain reliable orderings of ground and excited states. 
We find that the most stable fee phase is a spiral with an equilibrium lattice constant 3.56 A. 

PACS numbers: 71.15.Ap, 71.15.Mb, 75.50.Bb 
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I. INTRODUCTION 

Advances in experimental setups along the past two 
decades have allowed to grow in a controlled way, and 
characterize, tiny structures and devices, paving the way 
for the slow development of those fields of Materials Sci- 
ence now covered under the umbrella of Nanoscience. 
Ever more powerful computers and canny algorithms 
have also allowed to simulate larger and larger clusters of 
atoms or molecules, filling partially the bridge between 
theory and experiments. 

Molecular Dynamics packages based on Density Func- 
tional Theory (DFT)^ represent a specially useful set of 
theoretical tools in the analysis of such materials and de- 
vices. Among them, those which use basis sets (BS) of 
localized wave functions are particularly attractive since, 
one the one hand, can be easily adapted to work with 
order-N algorithms, and on the other, they can be writ- 
ten in a tight-binding language, which allows for easier 
analytical approaches. 

SIESTA is a simulation package which literally im- 
plements the tight-binding philosophy^ Indeed, it uses 
Norm-Conserving pseudopotentiala^. to integrate away 
core energy levels and very flexible BS made up of nu- 
merical atomic-like wave-functions to handle valence elec- 
trons. The use of minimal bases allow for fast calcula- 
tions which already provide a qualitative understanding 
of the simulated material. Accurate calculations can be 
performed at a higher numerical cost, using more com- 
plete BS. Assessment of the degree of reliability of a 
BS might be essential, since competing would-be ground 
states may in some instances have small energy differ- 
ences. Such analyses have already been performed for 
selected molecules and solids^, but not for magnetic ma- 
terials. Those studies show how both the number of wave 
functions used as well as their extent are variational pa- 
rameters, providing therefore a path for systematic im- 
provements of the accuracy of a simulation. A similar 



study for magnetic elements or materials seems to be 
highly desirable, since they have their own peculiarities 
and, in particular, are usually tougher to simulate. We 
have performed an exhaustive study of the degree of ac- 
curacy of the basis for Iron in most of its bulk phases 
as well as for small clusters. We find that SIESTA pro- 
vides a highly accurate description of the systems we have 
scrutinized, provided that a large numbers of extended 
orbitals is used. We will show below how BS regarded 
as fairly complete for simpler elements can provide dis- 
astrous results for this transition metal. 

The SIESTA package has built into it the possibility 
to cope with non-coUinear commensurate spin structures 
but only in the Local Density Approximation (LDA) . We 
have therefore coded the ability to compute non-coUinear 
arrangements of the spin moments in the Generalized 
Gradient Approximation (GGA), since LDA fails to pro- 
vide adequate ground states and lattice constants of a 
number of magnetic transition metals. Moreover, we 
have included the possibility to simulate non commen- 
surate spiral structureaSiS. 

While bcc (a— )Iron is firmly established to be a ferro- 
magnet, the scenario for fee (7— )Iron is more complex, 
since it stands at a crossing point between high spin ferro- 
magnetic and antiferromagnetic states, and the actual re- 
alization depends sensitively on its actual atomic volume 
and, possibly, strainsiS. Tsunoda discovered a decade 
ago that 7— Iron could be stabilized as pellets of radii up 
to 100 nm, with a lattice constant of 3.577 Aii-iS. He 
also found that the magnetic structure of the pellets was 
helicoidal, with pitch vector q^xp = (0.12,0, 1). A num- 
ber of authors have subsequently looked for theoretical 
low-energy coUinear and non-coUinear states appearing 

in such 7-phase.6,13il4,15,16,17,18,19,20_ 

Knopfle and co-workersS<, who used GGA or LDA and a 
full-potential implementation of Density Functional The- 
ory found that the ground state was indeed a spiral 
with the correct pitch vector. But since the Augmented 
Spherical Wave (ASW) method tends to overestimate the 
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atomic volume, it is difficult to extract what their equi- 
librium lattice constant oq might be. 

We have found a spiral ground state with gq of 3.56 
A, in excellent agreement with the experimental data of 
Tsunoda. We find two local minima when we plot the 
energy of the spiral state as a function of the pitch vector, 
at qi = (0,0,0.6) and q2 = (0.12,0,1), as other authors 
didii. We label these two states as SI and S2. The state 
SI is the global minimum for lattice constants down to 
3.47 A, and only below it is S2 the ground state. 

We believe that our results represent a significant 
methodological advance for the simulations of magnetic 
systems using the SIESTA package since, on the one 
hand, we lay down a firm ground for the reliability of 
atomic bases of different sizes and, on the other, we al- 
low for the description of interesting non-coUinear and 
spiral structures of the spin. 

The layout of this article is simple. Section II provides 
the theoretical backbone of the article; Subsection Ila is 
a brief reminder of the non-coUinear formalism as is ap- 
plied to DFT; we present subsequently details of our im- 
plementation of non-coUinearity and non commensurate 
spiral arrangements of spin for such structures. Section 
III is devoted to show and discuss our results for the sta- 
bility of the different states of Iron, both within LDA and 
GGA, using different BS, up to an optimal choice. We 
finish the article with a short summary. 



II. THEORETICAL BACKBONE 
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3. The Total Energy is a functional of the density 
matrix, E[h{r)] = T[h{r)] + EH[n{r)] + E^dHr)] 
which, upon variation provides with the effective 
single-particle Hamiltonian 



Hdft = ( -T^V^ + VH[n{r)] ] I + v,,[fiir)]. (4) 
Zm 



T and Eh are the kinetic and Hartree energy func- 
tionals, while 



E.^^ = j drfUn.m,u„.,Wn^Wm^^u^) (5) 

takes account of exchange and correlation, vh and 
Vxc, in Eq. (4), are the corresponding potentials. 

4. The spinor eigenfunctions of iJ, i^iir) can be 
used to compute the density matrix, since fi{r) = 
Tlii V'i "04 Each eigenfunction can individually be 
rotated to bring it back to coUinear form 



We devote this section to provide details of our im- 
plementation of non-coUinear GGA and of unconmcnsu- 
rate spiral arrangements. We believe it is useful to sup- 
ply first a backbone of non-coUinear DFTi&, which will 
help us discuss similarities and differences with the latest 
approaches^. 
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5. The exchange and correlation potential matrix, 
which is obtained by functional differentiation 
of the exchange and correlation energy, can be 
uniquely decomposed in terms of Pauli matrices, 



A. Brief presentation of non-collinearity in DFT 

1. In a non-coUincar material, the direction of the 
magnetization vector rh(f) — m ■ Um changes at 
each place in the sample, according to the angles 6*, 



Um = ( sin{9) cos{(f)), sin{9) sin{(j))^ cos{9) ) (1) 

2. The density matrix can be decomposed in terms of 
the electronic density n{r) and m(r) as 



n{r) = — (n{r) + 'm{r) ■ f) 



(2) 



where f denote the three Pauli matrices. There 
is a single rotation matrix, U{9{f),(f>{f)) = 
giTye/2 ^iT^ip/2 .^^jjich brings h{r) to collinear form. 
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where Vg = tr{vxcl)/'^ and Vy — tT{vxcT)/2. 



6. In LDA approximation, Vy is a function of only one 
vector, Um so that it must be proportional to it. It 
is then easily shown that Vxc is diagonalized by the 
same rotation matrix as ri. 



^^ ' ux 



U 



(8) 



where v^ = v„ [n^] . This can be interpreted phys- 
ically as rotating the whole system into a collinear 
reference frame, where Va{r) can be computed as 
in conventional LDA. 
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FIG. 1: Total Free Energy as a function of (a) number of k 
points in half the Brillouin zone (with a Grid cutoff of 400 Ry) 
and (b) Grid cutoff (with 4.000 k points in half the Brillouin 
zone). An optimized Double-Zeta (DZ) basis was used. 
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where A = (aa:,aj,,0), both approaches are analytically 
equivalent. Kleinman and BylandeB^SiS, added a spin 
stiffness term to their LDA exchange and correlation po- 
tential from which they obtained a spiral ground state 
with q— (1,0,0) at the atomic volume of bulk copper. 

The effective LDA Hamiltonian does not commute 
with the Pauli matrices fx,y unless the system be para- 
magnetic or collinear {0 be equal to 0°). Therefore the 
expectation values of Sx,y are not conserved in the itera- 
tive self consistency process of DFT. States with = 90" 
can be shown to be metastable, and therefore is con- 
served in this case. 



7. The GGA expression for the exchange and correla- 
tion energy contains also the vector V'^Um- There- 
fore the GGA potential includes both spin stiffness 
and antisymmetric exchange terms (Dzyaloshinskii- 
Moriya^iiSa) 



= Vsi+iVmUm+Vgrad'^'^Um+VcrossUmX^'^Um)f (9) 

This implies that v^c can not be fully diagonalized 
by the U rotation matrices. 



B. Description of Non-collinear commensurate and 
spiral states for a localized BS 

A convenient variational wave function for either 
molecules (or solids) with non-coUincar (commensurate) 
magnetic moments is 
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(14) 



UVxcU'' = vTJ'' + U{Vgrad^'^U^ + VcrossUm X V^U,„)TC/'f 

(10) 

An accurate description of non-coUinearity must in- 
clude a spin stiffness terrain, as is the case for classical 
localized spins. For instance, the Heisenberg Hamilto- 
nian is 



For solids, the above wave function can be easily 
rewritten so that it explicitly satisfies Block theorem. 
For helicoidal arrangement of spins of pitch vector q, the 
DFT Hamiltonian commutes with the operator T(i?, q) = 
U{0,q- R)T{R), which translates by a lattice vector R 
and then rotates about the z-axis. Since a wave function 
of the kind 
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where p is the spin stiffness. 

Our recipe to partly take account of gradient correc- 
tions is to neglect the stiffness contribution and evalu- 
ate only "5^°". We rotate the density matrix to bring 
it into collinear form, n'^°". We then compute its gra- 
dients Vn'^"" in this collinear reference frame, and the 
collinear potential, w^"" . We finally rotate back. Knopfle 
and coworkers chose an apparently different approach^. 
They rotated both h{r) and its gradient, Vn(r), and then 
evaluated v^c with the diagonal terms of both matrices, 
discarding the non-diagonal terms of Vn(r). But since 
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is an eigenfunction of T{R, q), a generalized Bloch theo- 
rem holds.* It must be stressed that such wave function 
is not an eigenstate of Hoft since, as noted above, a ro- 
tation by a constant angle 6*0 about the y-axis, U{6q,Q), 
does not commute with the Hamiltonian, unless ^o = 0- 
We have checked numerically that for any wave-function 
of the form above, the angle 9o is indeed not conserved 
by the application of Hdft, unless 9 = 90*^, which cor- 
responds to a metastable situation. Such wave function 
must therefore be regarded as purely variational. 
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FIG. 2: Energy difference (a) and number of iterations (b) 
as a function of the electronic temperature. The calculation 
was performed with 4.000 k points, a grid cutoff of 400 Ry 
(27.000 points) and an optimized DZ basis. 
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FIG. 3: (a) Evolution of the Free Energy of the three most 
stable states of Iron as a function of the size of the BS. AF 
= antiferromagnetic, HS = high-spin, FM = ferromagnetic, 
(b) Cohesive energy curves of those same three states, for the 
two minima found using the BSD. The discontinuous curves 
correspond to radii of 6 a. u. and the continuous to 10 a. 
u.The calculations have been performed using the GGA ap- 
proximation. 



III. RESULTS 

A. Preliminaries: choice of pseudopotential and 
integration grids 

SIESTA uses norm-conserving pseudopotentials 
(NCPS) « optimized so that their local part be smooth^. 
Izquierdo and coworker^4 proposed to generate the 
NCPS for Iron from the atomic configuration [Ar]3d^4s^, 
with core radii for 4s, 4p, 3d and 4f orbitals set equal 
to 2.00 a.u. They found that the optimal radius for 
partial-core corrections was 0.7 a.u. We have looked 
carefully for a pseudopotential which could produce a 
better fit to bulk bcc Iron. Our first criterium to de- 
termine the cutoff radii was to compare the eigenvalues 
of the valence shells of atomic Iron obtained from the 
NCPS and the all-electrons Hamiltonians and try to 
minimize their differences. We found that the radii 
obtained using such procedure were very different from 
each other. Moreover, they produced fits to the bulk bcc 
phase of poor quality as compared with the proposal of 
Izquierdo et al. Inclusion of the partial-core 3p levels 
into the valence did not help, mostly due to the fact 
that 3s electrons are still taken as part of the core and 
therefore there remains a strong overlap between valence 
and pseudo-core charge. 

SIESTA performs Brillouin zone integrations on a grid 
of Monkhorst-Pack special points typically extended to 
cover half of 1*2*2^. The integrand is also smeared by a 
Fermi function. Hamiltonian matrix elements are partly 
computed on a real space grid, whose fineness Ax is con- 



trolled by a grid cutoff, E^ « {tt/Ax)^/2. 

Since the energy of the different states need not shift 
rigidly when increasing accuracy and moreover compet- 
ing ground states for (fee) Iron have energy differences as 
tiny as 5 meV, we decided to set the number of k points, 
the electronic temperature and the grid cutoff to match 
an accuracy of about 1 meV. Figs. ^ and HI show typical 
results for the convergence of the energy of bcc ferro- 
magnetic Iron as a function of those parameters. All the 
calculations shown there were performed using the GGA 
functional as parametrized by Perdew and coworkers^^ 
and optimized either Double-Zeta or Triple-Zeta bases 
(see below). We find then that we need 4.000 k points 
and up to 700 Ry (which corresponds to 50.000 points in 
the real space grid) to meet the desired accuracy. Fig. 
121 shows that increasing the temperature to an optimal 
value of about 200 or 300 K speeds the convergence of 
the selfconsistent process significantly while not damag- 
ing the accuracy required for the energy. 



B. Optimization of the atomic basis and Phase 
diagram for bulk Iron 

SIESTA allows for a great flexibility in the use of BS of 
wave functions which describe valence electrons. For each 
species of atom, one may specify one shell of s, p, d and 
f orbitals. Within each shell, one may choose how many 
wave functions having the required angular symmetry are 
needed. A Single-Zeta basis (SZ) is equivalent to choos- 



ing just one. Completion of the basis leads to Double- 
Zeta and Triple-Zeta bases (DZ, TZ). In addition, one 
may polarize an orbital (P), which means adding wave 
functions which correspond to one higher angular mo- 
mentum unit^'. The minimum basis required to accom- 
modate the eight valence electrons of Iron would be SZ 
for both s and d orbitals, which provides a total of 6 wave 
functions per spin. SIESTA is set up to the default max- 
imum BS TZTP, which corresponds to 72 wave functions 
(WF) per spin. Using more wave is equivalent to filling 
up the Hilbert space and provides a better variational 
estimate of the ground state. Junquera and coworkers^ 
pointed out that the confinement radius of each orbital 
are also variational parameters. Very fast calculations or 
simulations of a large number of atoms may therefore be 
performed by the use of a Single-Zeta basis of rather con- 
fined orbitals. Such calculations usually provide much of 
the features of a material or device. But they are usu- 
ally regarded are pretty inaccurate, and DZ bases with 
polarized s orbitals are rather used. 

We have minimized a few BS ranging from SZ-SZ-SZ (9 
WF) to TZTP-TZTP-TZTP (72 WF). Fig. (a) shows 
how the convergence of the energy for the three most 
stable states of bulk Iron, e. g.: bcc ferromagnetic, fee 
ferromagnetic high spin and fee antiferromagnetic as a 
function of the number of orbitals used in the BS. We 
find that a TZ-TZ-TZ BS (27 WF) is essentially con- 
verged for p, d and f orbitals, since the Free Energy of 
the three states changes only a little if we polarize this 
BS. While we have not checked explicitly that a fourth 
Zeta for the s-orbital may still change somewhat the en- 
ergy, a inspection of the curve induces us to believe that 
our results are completely converged. It is also appar- 
ent that the states do not shift rigidly upon improving 
the accuracy. We shall see below that such an effect is 
specially damaging in the LDA approximation. 

We have paid special attention to the minimization of 
the BS DZP-SZ-DZ (BSD) and TZ-TZ-TZ (BST), where 
we have used a grid software program to look for local 
minima of the energy as a function of the radii of the first 
Zeta of s, p and d orbitals. We find a first local minimum 
for somewhat confined radii of about 6 a.u. and a deeper 
one for radii of about 10 a.u. We found that the energy 
still decreased upon looking further away, but thought it 
worthless to attempt to look for such next minimum. Fig. 
13 (b) shows that extended radii improve both the energy 
and the lattice constant substantially. For instance, the 
lattice constant of the bcc ferromagnetic state obtained 
using BSD in GGA approximation, as predicted by the 
first minimum is 2.90 A, while the second one gives oo = 
2.88 A. 

We find that DZP BS predict erroneous orderings of 
the ground and first excited states. Such effect is partic- 
ularly apparent in the LDA approximation. Fig. ^ (a) 
shows that BSD erroneously predicts that the LDA ferro- 
magnetic bcc state is more stable that the paramagnetic 
fee one. Fig. ^ (b) shows that usage of more complete 
BS correct such a mistake. Under such proviso, SIESTA 
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FIG. 4: Free energy and magnetic moment of the ground 
and lowest excited states of bulk Iron as predicted by LDA 
approximation, using (a) BSD and (b) BST, as a function of 
the atomic volume. 



provides pretty accurate results for the LDA predictions 
of the different physical magnitudes. For instance, the 
lattice constant, magnetic moment and bulk modulus of 
ferromagnetic bcc Iron are found to be 2.76 A, 2. OS/is 
and 2.68 Mbar, which compare extremely well with the 
best all-electrons Plane-wave calculations «. Moreover, 
the lattice constant for paramagnetic fee, 3.38 A, is also 
very similar to all-electrons estimate of 3.375 A, while 
the energy difference between both states is somewhat 
underestimated (55 versus about 70-80 meV)iSl. 

DZP Basis sets also provide awkward results in the 
GGA approximation, even though the relative stability of 
the lowest energy states is correct now, see Fig. El Never- 
theless, we find that the shape of the energy curves of the 
fee states change significantly when we increase the size 
of the basis from BSD to BST. Now, since the spiral state 
smoothly interpolates between the ferromagnetic HS and 
the antiferromagnetic ones, we have increased further the 
size of the BS. We have included more polarization or- 
bitals of p, d and f symmetry, and have found that the en- 
ergy of the three curves is essentially converged (see Fig. 
OJ. We find equilibrium lattice constant, magnetic mo- 
ment and bulk modulus of 2.85 A, 2.31fj,B and 1.83 Mbar 
for the ground state, which compare reasonably well with 
former all-electrons or ultrasoft-pseudopotentials-based 
plane waves calculationsSii^. We have computed the 
properties of the spiral state for lattice constants well 
below 3.54 A, so that the binding energy curve is a clear 
parabola, with a minimum at a = 3.56 A, very close to 
the experimental value (2.577 A). 

We have also simulated clusters with a number of 
atoms ranging from 2 to 5, using a BST and non 
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FIG. 6: Free energy and magnetic moment of the ground 
and lowest excited states of bulk Iron as predicted by LDA 
approximation, using (a) BSD and (b) BST, as a function of 
the atomic volume. 



FIG. 5: (a) Free energy and (b) magnetic moment as a func- 
tion of pitch vector q of the spiral state for different lattice 
constants ranging from 3.48 Ato 3.61 A. 



TABLE I: Bond lengths a (A), binding energy per atom E^ 
(eV/atom) and total magnetic moment M (/is) for Iron clus- 
ters up to 5 atoms calculated with a Triple-Zeta basis and 
GGA. 
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coUinear GGA, as shown in Table Q] Our calcu- 
lations compare very well with previous theoretical 
simulations^^iSSiSSiSSiSi and even improve them when 
comparisons are made with the experimental values 
found for the Fe2 cluster—. 



curves for lattice constants equal or larger than 3.58 A 
only have the qi minimum. The second minimum ap- 
pears when we decrease a at or below 3.54 A, becoming 
lowest in energy at a « 3.50 A. The curves correspond- 
ing to smallest lattice constants are very shallow. Their 
two minima have almost the same energy up to 1 meV, 
and are separated by energy barriers as small as 4 meV. 
If LDA were accurate enough for Iron, one would expect 
both phases not only to coexist but also to change dy- 
namically from one to the other. 

We finally discuss our results for the spiral structures in 
GGA, where we find that the S2 state has already clearly 
developed when a — 3.52 A, but that the ground state 
is SI down to lattice constants of 3.47 A. It can be seen 
again that the energy curves change rather much when we 
increase the size of the basis. One of the reasons is that 
the shape of the binding energy parabola also changes 
significantly (see Fig. EJ. But for BST, the ferromagnetic 
state has always considerably higher energy and there is a 
clearer asymmetric double- well structure with activation 
barriers of about 5-7 meV. 



C. Spiral states in the 7 phase 

We turn now to the predictions for the spiral state 
in LDA. We scan the energy as function of pitch vector 
along the FX and XW directions, where we find the two 
minima tfi and 92 we talked about in the introduction. 
On closer inspection of Fig. |S1 we see that the energy 



Marsman and Hafner have also performed simulations 
of 7-lron under tctrahedric, orthorombic and monoclinic 
distortions. However they obtained for the undistorted 
case a spiral state S2 with equilibrium lattice constant 
a = 3.49 A, much smaller than the experimental one. 
They also found that the equilibrium lattice constant for 
SI was a = 3.51 A. On the contrary, we obtain a state 
SI with a = 3.56 A, much closer to experiments. 
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the Iron phases, since for LDA it predicts that the bcc 
phase is less stable than the foe, which is contrary to all 
previous calculations. However, when we use a triple- 
zeta basis the results change dramatically for both LDA 
and GGA. We obtain for the ferromagnetic bcc and equi- 
librium lattice constant of 2.85 A, a magnetic moment of 
2.31 fiB and a bulk modulus of 1.83 Mbar, in excellent 
agreement with experiments. We have also simulated 
Iron clusters and we found a better estimate of the prop- 
erties of these materials than previous works. 

Finally, we have also made a profound study of the 7 
phase and we found an equilibrium lattice constant of 
about 3.56 A, closer to the experimental value of 3.577 
A, but with a spiral state SI instead of the experimental 
S2. However these simulations agree with previous works 
and even improve them. 



FIG. 7: (a) Free energy and (b) magnetic moment as a func- 
tion of pitch vector q of the spiral state for different lattice 
constants ranging from 3.48 Ato 3.61 A. 



IV. SUMMARY 

We have performed a thorough study of the different 
Iron phases in order to provide a good basis for future 
calculations. We have paid special attention to the min- 
imization of the different parameters used in the LCAO- 
DFT SIESTA code. We found that a grid cutoff of 700 
Ry, 4.000 k points and temperature smaller than 300 K 
are needed to meet an accuracy of about 1 meV. 

We showed that a double-zeta basis is not accurate 
enough to supply the correct ground state ordering of 
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